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Abstract 

From the non-equilibrium critical relaxation study of the two-dimensional 
Ising model, the dynamical critical exponent z is estimated to be 2.165 ± 
0.010 for this model. The relaxation in the ordered phase of this model 
is consistent with exp(—^t/r) behavior. The interface energy of the three- 
dimensional Ising model is studied and the critical exponent of the correlation 
length v and the critical amplitude of the surface tension oq are estimated 
to be 0.6250 ± 0.025 and 1.42 ± 0.04, respectively. A dynamic Monte Carlo 
renormalization group method is applied to the equilibrium properties of the 
three-dimensional Ising model successfully. 
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1 Introduction 



The study of the non-equilibrium relaxation turned out to be useful to study 
the critical-slowing down and other dynamical aspect of the Ising models JL), 

IJ,1,S1- 

For the critical relaxation, it is shown that the analysis based on the dy- 
namical finite-size scaling theory [[F] produces precise estimation of the value 
of f3/zv, where j3 and v denote the equilibrium critical exponents of the 
spontaneous magnetization and the correlation length, respectively, and z 
denotes the dynamical critical exponent. It is estimated directly from the 
initial relaxation curve of the magnetization m(t) which follows asymptoti- 
cally t~~^l zv . This method is simpler than the other previous methods which 
require the calculations of the two-time correlations and the estimations of 
the correlation times from them. In this paper, the values of z for the two- 
dimensional Ising model is studied from the non-equilibrium relaxation. The 
exact values of (3 and v of this model are known to be 1/8 and 1, respectively, 
and the value of z is estimated from the value of (3/zv. 

The non-equilibrium relaxation analysis is useful also for the relaxation 
process in the ordered phase which is not yet fully understood. It has been 
pointed out that the equilibrium relaxation of the magnetization may be 
stretched exponential, that is, 

C m {t)~e~ at \ (1) 

where C m (t) denotes the equilibrium time-correlation function of the mag- 
netization and a < 1. Two phenomenologies predicts a = 1/20 and 1/30 
for the two-dimensional Ising model. In this paper, the non-equilibrium 
relaxation in the ordered phase is studied assuming the similar stretched 
exponential decay, that is, 

m(t) -m s ~ ae~ bta . (2) 

The value of (3/zv for the three-dimensional Ising model was estimated 
to be 0.250 ± 0.002 J5|. The value of (3/v of this model has been estimated 
by several methods and the non-equilibrium critical relaxation provides one 
method to study its value. When the scale transformation with scale I is 
applied n times repetitively, the scale transformed magnetization is expected 
to follow the scaling form 

m{t,n) = l n ?/ u f{tl- nz ) (3) 
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at the critical point for large tl 



In this scaling region, the magnetization at the fixed time follows [10. I I 



m 



r p/v . (4) 



When we use the total magnetization M instead of the magnetization per 
site, it behaves as 

This relation provides a new method of Monte Carlo renormalization group. 
The study of /3/u for three-dimensional Ising model based on this scaling 
argument is also included. 

The second problem of this paper, the interface energy and the surface 



tension []12[. In the recent studies of this problem flT3|, |H| , |T5| , [376] , |T7| , |I8| , it 
seems to be believed that the standard Monte Carlo simulation with single 
spin updating dynamics is not an appropriate method for this problem. In 
this paper, however, it is shown that the naive simulation works successfully 
and a precise study of the interface energy is made. We reconfirm the values 
of the critical exponent and the critical amplitude, which is important to 



verify the universality hypothesis about the critical amplitude |L9]. 

The non-equilibrium relaxation study of the two-dimensional Ising model 
is described in the next section and the interface energy of the three-dimensional 
Ising model in the third section. The /3/u analysis is given in the appendix A. 
The last section contains the summary. 

2 Two-dimensional Ising model 

The non-equilibrium relaxation phenomena is studied for the two-dimensional 
Ising model in this section. 

2.1 Dynamics and algorithm 

The Hamiltonian of the ferromagnetic Ising model on square lattices is 

-(3H({a}) = K ]T a t a 3 (K > 0, a t = ±1), (6) 

\i-j\=l 



3 



where (3 = 1/kT and the summation runs over all the nearest neighbor site 
pairs. The dimensionless constant K is used as the inverse temperature in 
the following. The Ising spins at the lattice points are denoted by 

a (i,j) (* = 1 , 2 ,-'-, L x,j = 1,2, (7) 

where L x and L y are assumed to be L + 1 and L, respectively, and L is 
an even integer. Skew boundary condition is applied to L x direction and 
periodic boundary condition to L y direction. There is no intrisic dynamics 
of this model and the dynamics studied here is the single-spin updating 
stochastic dynamics with two-sublattice flip sequence and Metropolis-type 
transition probability. The spins are classified into two sublattices defined 
by 

VLa = i + j = even } and Qb = {^(ij)] i + 3 — odd }. (8) 

The spins in VLa are updated firstly and then those in fl B are updated. The 
transition probability is 

PmOo -> -o- ; (Ji, a 2 , a 3 , a 4 , ) = min{l, exp[-2 J E(a ; a u a 2 , a 3 , a 4 )]}, (9) 

where <Ti, <j 2 , <t 3 and cr 4 are the nearest-neighbor spins of cr and the local 
energy E denotes 

E(a ; oi, a 2 , a 3 , cr 4 ) = Ka^ax + o 2 + (7 3 + a 4 ). (10) 

This dynamics is efficiently vectorized and it is the reason why it is used. 
64 independent systems with the same lattice are simulated simultaneously. 
The 64 spins on the same site of these 64 systems are coded in the same eight- 
byte integer variables, that is, the independent-system coding technique is 
employed |20], ^1], ^2|. The spin-updating procedure is realized by 15 logical 



operations and the generation of one random integer. The condition that a 
spin cr is flipped to — a is 

r <Pm(<7o,-<7o;<7i,<72,<7 3 ,<74), ( U ) 

where r is a random number. The transition probability pm is the function of 
local energy E(a Q ; <Ji, <r 2 , &3, cr 4)- The spin variable <7j which takes —1 or +1 
value is stored in one bit of one computer word which takes or 1 value and 
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the corresponding bit is denoted by I{. Now the above condition is expressed 
by 



h ® Io + h ® I Q + h ® h + h ® h + 2 * IXlT(r) + IX2T(r) > 2, (12) 

where © denotes the exclusive-OR operation and IXlT(r) and IX2T(r) are 
the functions of random number r defined by 

IXlT(r) { q otherwise 

and 

1 e -8 ^ < r < e _4if 
IX2T(r) { q otherwise 

This condition should be checked by the bitwise logical operations. The 
following operations realize it with 14 logical operations. In the following, 
< J, J > denotes 2* I + J. 

1. Make Jl and J2 so that < Jl, J2 >= II + 72 + 73. It takes five 
operations. 

2. Take the exclusive-OR of Jl, J2 and J4 with 10, respectively. It takes 
three operations. The same notations are used for the results. 

3. Make Kl and K2 so that < Kl, K2 >= 14 + IXlT(r) + IX2T(r). It 
takes three operations. 

4. Make ID which holds the truth value of spin-flip condition by 

Kl or Jl or (J2 and K2). (15) 
It takes three operations. 

Further details of this updating algorithm are found in the FORTRAN rou- 
tines given in the Appendix B. To do the independent simulation of the 
systems at the same temperature with the independent-system coding tech- 
nique, the recycle algorithm of the random numbers is used by shuffling the 
Boltzmann factor tables|2^]. The RNDTIK routine|24| are used for the ran- 
dom number generations. 
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The NEC SX3/11 of Koln university and the MONTE-4 of Japan Atomic 
Energy Research Institute are used for the simulation. The performance of 
the simulation codes developed for the present study, named IS2DVP, on 
these vector processors is shown in Table [1]. The MONTE-4 has four pro- 
cessors and they share the main memory. So the simulation speed on one 
processor is influenced by other jobs on other processors. The performance 
in Table [l] is measured with/without other jobs on other processors. The 
concurrent jobs were Monte Carlo simulations of two- and three-dimensional 
Ising models which are highly vectorized. The speed-down shown in Table [l] 
under multi-job environment can be regarded as the worst case. If the concur- 
rent jobs are badly vectorized jobs or the jobs which are highly vectorized but 
do not generate the main memory access massively, the speed-down becomes 
smaller. In fact, almost no speed-down was observed when the performance 
of IS2DVP is measured with three lowly vectorized jobs. 

The value of magnetization at tMCS from all up initial condition is de- 
noted by m(t). This time-dependent magnetization of the above-mentioned 
dynamics is estimated by the simulation. It is estimated from the values of 
magnetization at t of independent Monte Carlo runs. 



2.2 At the critical point 

The critical point of the square-lattice Ising model is K c = J/ k B T c = 0.440686793 • • • 
where the time-dependent magnetization m(t) decays to zero following asymp- 
totically to 

m(t) ~ t~ p/zv = r l/Sz , (16) 

where the exact values of the critical exponents of the spontaneous magneti- 
zation and the correlation length, (3 = 1/8 and v = 1, are used. 

The simulations are listed in Table ^. The deviation of m(t) due to the 
finite lattice is expected to be exponentially small for large system. The 
system-size dependence at several time are shown in Fig. [I]. No size depen- 
dence is observed up to t = 400MCS, although it is observed at t = 500MCS. 
In the following, the estimated values of m(t) for L = 1500 lattice up to 
400MCS which are shown in Fig, ^ are used for estimating the value of z. 

Firstly the estimator 

R(t,m)=t[ m ^~ 1) -1] (17) 
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is used for this purpose. The average of this R(t) over N successive t, that 
is, 

i t+N-l 

R%ve.(r(t,N),m) = - £ R(r,m), (18) 

iV T=t 

is used to reduce the statistical fluctuation, where the mean time, r(t, N) = 
t + (N - l)/2. The values of zf ve _ = l/(8Rf ve .) is plotted in Fig. |. If 
the form of m{t) is t~ x (a + a\jt + a 2 /t 2 + • • •), -R^ve. ( r > m if)) is the form 
of A + di/r + d 2 /r 2 + • • •. It is observed that the value of z converges to 
2.165 ± 0.005 assuming that the correction to the asymptotic simple-power 
law is analytic, which is consistent with our data. 

The estimates for z using the least-square fitting analysis are shown in 
Fig. |j. The function of the form of logm(t) = a — blogt is fitted to the 
successive fifty values of magnetization, that is, m(t), m(t + 1), ■ ■ • m(t + 48) 
and m(t + 49). The values of zj^j n g = 8/6 are plotted. In the large time 
limit, the value of z estimated from this Fig. |4] is consistent with the estimated 
z = 2.165 ± 0.005 from the ratio analysis. 

There are many theoretical studies on this value of zpSL |26|, E7], 



The estimates for the two-dimensional Ising model are listed in Table y 
and they are plotted in Fig. ^|. It is observed that the simulational values 
which includes the estimates obtained by temperature dependence analysis, 
finite-size scaling analysis, Monte Carlo renormalization group analysis, non- 
equilibrium relaxation analysis and other analysis of the simulational results 
are converging. 

Above analysis assumes the form of the function m(t) to be a simple 
power-law with simple correction terms and the results might have unknown 
systematic error caused by this assumption. Therefore the naive error region 
might be underestimated and an overestimated error region will be safer as 
our present conclusion. So we adopt z = 2.165 ± 0.010 as our conclusion. 

There is a room for the qualitatively different assumption for the critical- 
slowing down phenomena from above-mentioned simple power law. For 
example, a behavior with logarithmic correction to the naive diffusive be- 
havior, z = 2, was proposed [|j5[] . Assuming m(t) to be t~ l3 ^ zu (\ogt)~ c , 
where /3/zu = 1/16, the values of t 1//16 m(t) is plotted in Fig. [5[ The value 
of c is roughtly estimated to be 0.04 from the tangent of the plot about 
t = 200 ~ 400(logt = 5.3 ~ 6.0), although the convex curvature seems to 
remain. 
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2.3 In the ordered phase 



The simulations listed in Table |] are made. Each run started from all up 
initial configuration. In the previous section, we observed that the value of 
the magnetization of L = 1500 lattice up to 400MCS at the critical point can 
be regarded as values in the thermodynamics limit with the present accuracy. 
Now the simulation is made at off-critical temperature in the ordered phase 
and the maximum Monte Carlo step is up to 90 or 35. So the results can also 
be regarded as the values in the thermodynamics limit within the present 
accuracy. 

The estimated values of magnetization are plotted in Fig. [7]. The exact 
solution of the spontaneous magnetization, 

m s = [l £ -] 1/8 , (19) 

L sinh 4 (2iT) J ' V ' 

is used. The asymptotic slopes of these curves correspond to the value of a. 

At t = where m — 1, the values of ln[— ln(m — m s )] are 0.647, 0.819, 
0.943 and 1.041 for K = 0.47, 0.49, 0.51 and 0.53, respectively, and these 
values correspond to hit = — oo. Therefore the tangents of the curves in 
Fig. [7| approach the final value from smaller side and the curves are concave. 

The value of a observed from the Fig. |7| is clearly larger than 1/3. The 
results for K = 0.51 and 0.53 are consistent with a = 1/2. For K = 0.47 and 
0.49, the curved do not reach a = 1/2 behavior yet, but the final tangent are 
already larger than 1/3 which is observed in Fig. ^. 

We cannot exclude the possibility that the curves continue to become 
steeper for longer time and finally the simple exponential decay(a = 1) ap- 
pears. We can conclude that the conjecture a = 1/3 is not true for the 
non-equilibrium relaxation of this model in the presently studied time re- 
gion. 



3 Three-dimensional Ising model 

The interface energy of the three-dimensional Ising model is studied in this 
section. 
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3.1 Monte Carlo simulation 



The interface energy is estimated by the difference of the energy in the sys- 
tems with and without a interface. In this study, the different boundary 
conditions are used to realize the constraint for the interface. One lattice 
has periodic boundary conditions for all of the x, y and z directions. In this 
case, the configuration without interface is dominant because the present 
simulation is started from all up initial configuration and we take lattice size 
to infinity at fixed t. The other lattice has anti-periodic boundary condition 
in the z direction and periodic to the x and y direction. This means that the 
lattice is pasted into torus with anti-ferromagnetic bond with coupling con- 
stant —K(< 0) for z direction. In this case, a single interface configuration 
is dominant because the all-up initial configuration has one interface in the 
case of the anti-periodic lattice. 

The lattice size is specified by two integers, L and H, in the following. We 
use the box of cubic lattice of the size of L x L x H for x, y and z directions. 
Odd L and even H are used for the purpose of vectorizing the Monte Carlo 
dynamics. 

The total energies of periodic and antiperiodic lattices are denoted by E p 
and Egp. The interface energy per unit area denoted by e s is defined by 

e s (K,L,H) = (E ap -E p )/L 2 . (20) 

In the thermodynamic limit, the critical behavior of this interface energy 
is expected to be 

e s (K) = e s (K, L = oo, H = oo) = ^(1 - ^f v ~\ (21) 

where Oq denotes the critical amplitude of the surface tension and the factor 
is added instead of 1/K C because it gives a better description of the 
present results and the critical behavior is essentially the same. 

The simulations listed in Table |5] were made on MONTE-4 of JAERI 
and NEC SX3/11 of Koln University. In each simulation 64 independent 
simulations were run simultaneously using the independent-system coding 



technique fl20|, pi] , p2| and the recycle algorithm for the random number |23 
Half of these 64 systems have the periodic boundary condition and the rest 
have the anti-periodic boundary condition. All the simulations here are made 
in the ordered phase above the surface roughening temperature. The triplet 
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(A, B, C) in the MCS column of Table denotes the length of the simula- 
tion and the energy measurement timing. From all-up initial configuration, 
A MCS are made without energy calculation. The value of A is selected 
to be twenty or more times larger than the correlation time of the energy. 
Then the energy is calculated C times in every B MCS and one estimate 
of the interface energy is calculated using the estimates of the energy of 64 
systems. These operations are repeated five times independently. If five run 
can be made within the CPU time limit of one batch job, the initialization of 
the configuration is not made for every run but the last configuration of the 
previous run is used as the next initial configuration. The obtained five esti- 
mates of the interface energy are used to estimate the average, e s (K, L, H), 
and its error. 

3.2 Finite size effect 

For a fixed L, the interface energy will be independent of the height H if 
H is larger than the interface thickness which is known to increase as logL, 
as far as there is no interface in the periodic lattice and there is only one 
interface in the anti-periodic lattice. Fig. ^] shows the H dependence of the 
estimates of e s (K, L, H). It is observed that e s (K, L, H) is independent of H 
if H is not very small. The horizontal solid line and the broken lines show 
the estimate of e s (K,L) = e s (K,L,H = oo) and its la error region, which 
is estimated from the data H < 20, that is, all the data except two small H 
points shown in the figure. 

For every L at every K, lattices with several H are made. The values of 
the interface energy were estimated from the estimates for the if-independent 
region. The region of H used for this procedure is given in Table |^ with the 
values of e s . 

Now the extrapolation to L — ► oo has to be made. The finite size cor- 
rection to € S (K) is expected to be the order of 1/L 2 from the capillary wave 
contribution, that is, 

e s (K,L) = e s (K,L = cx>) + ^, (22) 

up to the first order correction. The values of the interface energy in the ther- 
modynamic limit e s (K) = e s (K,L = oo) are estimated from the estimates 
of e s (K,L) with this size dependence. This 1/L 2 correction is appropriate 
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to our results. One example of this extrapolation is shown in Fig. [II]. The 
estimated values of e s (K) are also given in Table ||. 



3.3 Critical behavior 



The critical exponent v and the critical amplitude Oq are estimated from our 
estimates of e s . The present estimates of e s are plotted in Fig. |TT . 



The value of the critical point K c has been estimated by several methods 
and usually is in the region 0.22165 ~ 0.22166|82| and this accuracy is enough 
precise for the present analysis. So K c = 0.221657|S3L B4] is used for the 
present analysis. 

The tangent of the consecutive two points in Fig. [IT] are shown in Fig. |T^| 
and the value of 2v — 1 at K = K c is estimated to be 0.250(5) from this 
figure. This means that the estimate of v is 0.6250 ± 0.0025. This values 
of v is consistent with the result of a recent Monte Carlo renormalization 
study[|2 



The previous estimates of v are also cited in Ref. |p2 



The value of a is estimated from Fig. 13 which shows the values of 



K K 



(23) 



from our estimates of e s (K) for several values of v. From this figure and above 
estimate of v, we estimate that the value of <To is 1.42 ± 0.04. The estimates 
in the previous works are shown in Table and are mostly confirmed by our 
more accurate data. 



4 Summary 



First the non-equilibrium relaxation of the two-dimensional ferromagnetic 
Ising model was studied. The critical relaxation exponent z was estimated 
to be 2.165 ± 0.010. This value is smaller than the recent other estimations 
for this value from the same non-equilibrium relaxation method by Miinkel 
et al[H and from the high-temperature expansion study by Damman and 



Reger|79[. But the upward behavior of the two-dimensional data shown in 
Fig. 1 of Ref. can be interpreted as a sample fluctuation and this might 
cause a overestimation. The coefficients of the high-temperature expansion 



given in Ref. [79| were reanalyzed by Adler[]8T| and an estimate consistent 
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with the present result was reported. The relaxation process in the ordered 
phase is stretched exponential. The exponent is larger than 1/30 and it is 
consistent with 1/2||. 

Secondly the critical exponent v and the critical amplitude <j$ are esti- 
mated to be 0.6250±0.0025 and 1.42±0.04, respectively. The fully vectorized 
simulation with single spin updating dynamics can efficiently study the in- 
terface. When we used only about 20 hours of NEC SX3/11, the data were 
as accurate as those of Ref. [[H] from 10 3 work station hours. 

Finally the new Monte Carlo renormalization group method produced 
reasonable estimate of (3/u with a quite small scale simulation. 



The total number of updated spins are 3.19 x 10 
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Appendix A 

A new Monte Carlo renormalization group method is applied to the esti- 
mation of /3/u value of the three-dimensional Ising model. The value of 
the magnetization of the scale transformed configuration with scale I = 2 
majority-rule real-space scale transformation is estimated. A 256 3 lattice is 
simulated up to 20000MCS from all-up initial configuraion with an appro- 
priate Monte Carlo dynamics. Totally 45 independent simulations are made 
and therefore the total number of the updated spin is 1.51 x 10 13 . 
Following the eq. (|J), the values of X(n, t) defined by 

%|) J, 6 [M(„-M)/M M ] (M) 

are studied. This A(n) is expected to be d — (3/u, where d is the lattice 
dimensionality, in the limit of t — > oo if the infinitely large lattice is studied. 
The behavior of the A(n, t) is shown in Fig. [T4|. From this figure, the value 
of A(2, t = oo) is estimated to be 2.4625(1). In table |8|, the estimated values 
of A(n) = \(n,t = oo) are shown for 256 3 and 1024 3 lattice. The 1024 3 



12 



lattice was simulated up to 5644 MCS for one sample and the number of the 
updated spin is 6.06 x 10 12 . 

When we extrapolate these estimates of A(n) are extrapolated to n — > oo 
with 

\( n ) = A + d~ nuJ , (25) 

the value of A behaves as shown in Fig. [15|. The estimates for the 256 3 lattice 
given in Table |8| are used for the extrapolation. From the present preliminary 
data, the value of uj is estimated to be 0.7(2) and the corresponding value 
of A is 2.492(8). this means the value of fi/v = 0.508(8). If the estimate 
u = 0.80 ~ 0.85f8| is assumed, A is estimated to be 2.487 ~ 2.490 and 
P/v = 0.513 ~ 0.510. If we assume u = 1, A and 0/u are 2.486 and 0.514, 
respectively. 



Appendix B 

Two routines of IS2DVP code which are specific to the two-dimensional 
Ising model are shown. BFL2DF generates the Boltzmann factor tables and 
SU2DSK updates the spin configuration. 

C IS2DVP: Monte Carlo simulation of two-dimensional Ising Model 
C 

C Boltzmann-f actor table preparation routine 
C METROPOLIS-TYPE TRANSITION PROBABILITY 
C 1992.09.19. VERSION 1.00 BY N0BUYASU IT0 

C 

C TK(64) R*8 : Give a inverse temperature of 64 systems 

C 1X1(0 :IRLST) , IX2(0:IRLST) 1*8 : The BFT is generated in them. 

C IRLST is maxinum number of the random integer. 

SUBROUTINE BFL2DF (TK , 1X1 , 1X2 , IRLST) 
PARAMETER ( I WIDTH=64) 
DIMENSION TK(IWIDTH) 
REAL* 8 TK 

DIMENSION 1X1 (0 : IRLST) , 1X2 (0 : IRLST) 
DIMENSION IKl(IWIDTH) ,IK2(IWIDTH) ,I2P(IWIDTH) 
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REAL* 8 TNORM 

C 

TNORM=DFLOAT ( IRLST) 
DO 10 I=1,IWIDTH 

IK1 (I) =IDINT (TNORM*DEXP (-8 . ODO*TK (I ) ) ) 
IK2 (I ) =IDINT (TNORM*DEXP (-4 . ODO*TK (I ) ) ) 
10 CONTINUE 

DO 30 I=1,IWIDTH 

I2P(I)=ISHFT(1, (IWIDTH-I)) 
30 CONTINUE 

DO 60 1=0, IRLST 
IX1(I)=0 
IX2(I)=0 
60 CONTINUE 
*VDIR NODEP 

DO 50 I=1,IWIDTH 
DO 40 J=0, IRLST 
IXT1=0 
IXT2=0 

IF(J.LE.IK1(I))IXT1=I2P(I) 

IF(J.LE.IK2(I))IXT2=I2P(I) 

I2=IE0R(IXT1,IXT2) 

I1=IAND(IXT1,IXT2) 

IX2(J)=I0R(IX2(J) ,12) 

IX1(J)=I0R(IX1(J) ,11) 
40 CONTINUE 
50 CONTINUE 
RETURN 
END 

C IS2DVP: Monte Carlo simulation of two-dimensional Ising Model 
C 

C Spin configuration updating routine for 

C ferromagnetic Ising model with Metropolis-type transition 
C probabilities 

C BFL2DF routine should be used for BFT preparation 
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C 1992.09.19. VERSION 1.00 BY NOBUYASU ITO 

C 

C ISTEP 1*8 : number of Monte Carlo sweeps 

C LI, L2 1*8: lattice size 

C IS((-L1+1) : (L1*(L2+1)) 1*8: spin configuration 

C IRD(L1*L2) 1*8: work area 

C IRLST 1*8: Maxinum number of random integer 

C 1X1 (0: IRLST) , 1X2(0: IRLST) 1*8: Boltzmann factor tables (BFT) 

SUBROUTINE SU2DSK (ISTEP , LI , L2 , IS , IRD , IRLST ,1X1, 1X2) 

DIMENSION IS((-L1+1) : (L1*(L2+1))) 

DIMENSION IRD(L1*L2) 

DIMENSION 1X1 (0 : IRLST) , 1X2 (0 : IRLST) 

NSYS=L1*L2 

LS=L1 

DO 10 IMCS=1, ISTEP 
CALL RND02I(NSYS,IRD) 
IFIRST=1 
*VDIR NODEP 

40 DO 20 I=-LS+1,0 
20 IS(I)=IS(I+NSYS) 
*VDIR NODEP 

DO 30 I=NSYS+1,NSYS+LS 
30 IS(I)=IS(I-NSYS) 
*VDIR NODEP 

DO 50 IJ=IFIRST,NSYS,2 
I0=IS(IJ) 
I1=IS(IJ+1) 
I2=IS(IJ-1) 
I3=IS(IJ+L1) 
K1=IE0R(I1,I2) 
K2=IAND(I1,I2) 
J2=IE0R(K1,I3) 
K3=IAND(K1,I3) 
J1=I0R(K2,K3) 
J1=IE0R(I0,J1) 
J2=IE0R(I0,J2) 
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I4=IS(IJ-L1) 

I4=IE0R(I4,I0) 

IRT=IRD(IJ) 

IX1T=IX1(IRT) 

IX2T=IX2(IRT) 

K2=IE0R(I4,IX2T) 

K1T=IAND(I4,IX2T) 

K1=I0R(K1T,IX1T) 

ID=I0R(J1,K1) 

K4=IAND(J2,K2) 

ID=I0R(ID,K4) 

IS(IJ)=IEOR(IS(IJ) ,ID) 
50 CONTINUE 

IF ( IFIRST . EQ . 1 ) THEN 

IFIRST=2 

GOTO 40 
END IF 
10 CONTINUE 
RETURN 
END 
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Figure 1: The size-dependences of m(t) at (a)t = 10, (b)100, (c)400 and 
(d)500 are shown. 

Figure 2: The estimated values of m(t) for L = 1500 lattice up to 400MCS 
are shown in logarithmic scale. They are regarded as those in the thermody- 
namics limit within the present accuracy. 

Figure 3: The value of the averaged local exponent ^|ve. ^ s shown. The 
magnetization obtained for L = 1500 lattice is used. 

Figure 4: The estimates of z from the least-square fitting are shown. r(t, 50) 
denotes the mean time of the data, that is, t + (50 — l)/2. 

Figure 5: The values of m(t) is plotted assuming the logarithmic correction. 

Figure 6: The estimates for z are plotted. The points marked with • and o 
show simulational and other results, respectively. The horizontal axis shows 
the publication year except the recent preprints. 

Figure 7: The non-equilibrium relaxation curves of magnetization are shown. 
The point marked by •, o, A and o correspond to K — 0.47, 0.49, 0.51 and 
0.53, respectively. The tangent of the sold line is a = 1/2. Those of broken 
lines are a = 1/3 (gentler one) and 1 (steeper one). 

Figure 8: The magnified figure of the magnetization at K=0.47. The solid 
line has the slope of 1/3. It is observed that the exponent a is already larger 
than 1/3. 

Figure 9: The height dependence of the interface energy is shown. The 
estimates of e s (K, L, H) for L — 21 are plotted as an example. 

Figure 10: The L dependence of the interface energy is shown. The estimates 
of e s (K, L) at K = 0.24 are plotted as an example. 
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Figure 11: The critical behavior of the interface energy is shown. Ke s (K) 
is plotted versus K — K c in the logarithmic scale. 

Figure 12: The values of the exponent 2v — 1 from the consecutive two 
points are shown. 

Figure 13: The values of <Jq{v) are plotted. *, A, •, o and o show the 
estimates for v = 0.6300, 0.6275, 0.6225 and 0.6200, respectively. 

Figure 14: The behavior of the A(2, i) is shown. Every point denotes an 
averaged value of A(2,t) over successive 1000 MCS. 

Figure 15: The extrapolated values of X(n) to n — > oo assuming the value 
of uj are plotted. 



machine 




performance 








U 


U+M 


one processor of 





1.31 


0.77 


MONTE-4 


1 


1.1 ~ 1.2 


0.6 ~ 0.7 




2 


0.9 ~ 1.0 


0.5 




3 


0.7 ~ 0.9 


0.4 ~ 0.5 


SX3/11 




1.05 


0.63 



Table 1: The performance of the simulation code, IS2DVP, is shown. The 
simulation speed in the case of 1501 x 1500 lattice is shown. The unit of the 
figures is GUPS, that is, giga(10 9 ) update per second. The figures in "U" 
and "U+M" denote the speed of the spin updating procesures only and of the 
spin updating procesure with magnetization calculation in every Monte Carlo 
sweep. The number A/j b in MONTE-4 shows the number of concurrent jobs 
on other processors of MONTE-4. The concurrent jobs are highly vectorized 
with dense memory access. Therefore the speed in 7Vj ob > 1 case shows the 
worst case. 
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L -^sample updated spins 

1000 1.54 x 10 5 7.69 x 10 13 

1500 1.18 x 10 6 1.33 x 10 15 

2000 6.14 x 10 4 1.23 x 10 14 

3000 1.54 x 10 4 6.91 x 10 13 



Table 2: The simulations made for the non-equilibrium critical relaxation 
study are listed, ^ sam pi e denotes the number of independent runs. The 
simulation up to 500MCS is made for every sample from all-up initial condi- 
tion. 
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year 


method 


1968 




1969 


HTE 


1971 


HTE 


1972 


RG 


1973 


MC 


1975 


RG 


1976 


MC 


1976 


HTE 


1976 


MCRG 


1978 


RG 


1978 


RG 


1978 


RG 


1978 


RG 


1979 


MC+FSS 


1979 


RG 


1979 


RG 


1979 


RG 


1979 


RG 


1979 


MCRG 


1980 


TM+FSS 


i you 


pn 

1 \ v r 


1980 


RG 


1980 


RG 


1980 


RG 


1981 


MC 


1981 


MCRG 


1981 


RG 


1982 


MC+NER 


1982 


MC+FSS 


1982 


MC 


1982 


MCRG 


1982 


MCRG 



z reference 

L75 !| 

2.00 ±0.05 g 

2 g 

2.0 g 
i.9o ±0.10 g 

1.82 g 

2.30 ±0.30 H 

2.125 ±0.01 H 

1.4 ±0.4 g 

2.2 g 

2.i9 gg 

1.7 g| 
1.67 H 

2.0 ±0.1 H 

2.09 @ 

1.96 |45| 

2.064(energy) g 
1.819(magnetization) |46| 

1.85 ±0.15 |47| 

1.99 m 

2.22 m 

i.8i9 g 

2.2 g 

i.79i g 

i.9 g 

2.22 ±0.13 || 

2.126 g 

2.1 g 
2.o g 

2.2 ±0.1 g 

2.23 g 

2.08 g 
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1 OfiQ 

lyoo 




9 1 9 


pl| 


lyoo 


DP 


9 9 






ruor,rl 1 Hi 


9 ^9 


E£l 


lyo4 




n 1,1 in no 

z.14 ± U.Uz 


J64| 


1 QS/1 




o 
1 


Jbo| 


lyoo 


DP 
lXXjt 


1. / O 


P b l 


lyoo 


DP 


l.oOO 


P' 


1 OS ^ 

lyoo 




9 9 


E°| 


lyoo 




9 1 q j- n nQ 
z.lo ± U.Uo 


P 9 I 






9 1 + n i 


IV 


1986 


MCRG 


2.19 + 0.05 


71 


1987 


MC+NER 


2.16 + 0.05 


72 


1987 


MC+FSS 


2.132 + 0.008 


73 


1987 


MC+FSS 


2.17 + 0.06 


74 


1988 


Q2R 


2.16 


75 


1988 


MC 


2.22 


75 


1988 


MC+NER 


2.076 + 0.005 


FBI 


1990 


HTE 


2.34 + 0.03 


77 


1992 


MC+NER 


2.2 


n 


1992 


HTE 


2.183 + 0.005 


79 


1993 


HTE 


2.21 + 0.005 




1993 


MC+NER 


2.21 + 0.03 


1 


1993 


HTE 


2.165 + 0.015 


Mil 



1993 MC+NER 2.165 ± 0.010 present study 



Table 3: The estimates of z are listed. The HTE, RG, MC, FSS and NER in 
this table are the abbreviations of high-temperature expansion, renormaliza- 
tion group, Monte Carlo, finite-size scaling and non-equilibrium relaxation, 
respectively. 
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K 


L 




■^sample 


updated spins 


0.47 


1500 


90 


1.02 x 10 6 


2.07 x 10 14 


0.49 


2000 


35 


3.58 x 10 5 


5.02 x 10 13 


0.51 


2000 


35 


3.58 x 10 5 


5.02 x 10 13 


0.53 


2000 


35 


1.54 x 10 5 


2.15 x 10 13 



Table 4: The simulations made for the relaxation study in the ordered phase. 
K and T max . show the inverse temperature and the Monte Carlo steps for 
each run. 
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L 


K 


H 


MCS 


11 


0.23, 0.24, 0.25 


6, 10, 16, 20, 26, 30 


(IK, 10, IK) 


15 


0.23, 0.24, 0.25 


6, 10, 16, 20, 26, 30, 36, 40, 46 


(IK, 10, IK) 


21 


0.23, 0.24, 0.25 

7 7 


10, 16, 20, 30, 40, 50, 60 


(2K, 10, 2K) 




0.23 


80, 100, 150, 200, 250, 300, 350, 400 


(2K, 10, 2K) 
(10K, 10 2 , IK) 




0.225, 0.2275 


50, 74, 100, 124, 150, 174, 200, 250 


25 


0.23, 0.24, 0.25 


10, 20, 30, 40, 50, 60, 70, 80 


(2K, 10, IK) 


31 


0.23, 0.24, 0.25 

7 7 


10, 20, 30, 40, 50, 60, 70, 80, 90 


(2K, 20, IK) 




0.235, 0.245 


50, 60, 70, 80, 90 


(2K, 10, 2K) 
(10K, 10 2 , IK) 




0.225, 0.2275 


50, 74, 100, 124, 150 


35 


0.23, 0.24, 0.25 


20, 40, 60, 80, 100, 120 

7 7 7 7 7 


(2K, 20, IK) 




0.235, 0.245 


50, 60, 70, 80, 90 


(2K, 10, 2K) 


41 


0.23, 0.24, 0.25 


20, 40, 60, 80, 100, 120 


(2K, 10, IK) 




0.235, 0.245 


50, 60, 70, 80, 90 


(2K, 10, 2K) 




0.225, 0.2275 


100, 124, 150 


(20K, 200, IK) 


45 


0.235, 0.245 


50, 60, 70, 80, 90 


(2K, 10, 2K) 


51 


0.23, 0.24, 0.25 

7 7 


26, 50, 76, 100, 126, 150, 176, 200 


(2K, 10, 2K) 




0.235, 0.245 


50, 60, 70, 80, 90 


(2K, 10, 2K) 




0.225, 0.2275 


100, 124, 150 


(20K, 200, IK) 


61 


0.23, 0.24, 0.25 


30, 60, 90, 120, 150, 180, 210 


(2K, 10, 2K) 




0.225, 0.2275 


100, 124, 150 


(20K, 200, 500) 




0.235, 0.245 


50, 60, 70, 80, 90 


(2K, 10, 2K) 


71 


0.225, 0.2275 


100, 124, 150 


(20K, 200, 500) 


81 


0.23, 0.24, 0.25 


40, 60, 80, 100, 120 


(2K, 10, 2K ) 




0.235, 0.245 


50, 60, 70, 80, 90 


(2K, 10, 2K) 


101 


0.225 


100 


(20K, 200, 400) 




0.225 


124 


(20K, 200, 300) 




0.225 


150 


(20K, 200, 250) 




0.23, 0.24, 0.25 


50, 60, 70, 80, 90, 100 


(2K, 10, 2K) 


151 


0.23, 0.24, 0.25 


50, 60, 70, 80, 90 


(2K, 10, 2K) 


171 


0.23 


50, 70, 90 


(2K, 10, 2K) 



Table 5: The simulations made for the interface study are listed. The 
three numbers in MCS column show the MCS schedules and their meaning 
is explained in the text (IK = 1000). Totally 1.26 x 10 15 spins are updated. 
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K 


L 


TT 

H 


A T 

N H 




o oo r 

0.225 


O 1 

31 


50-150 


5 


o rrnn/on\ 

2.5592(32) 




A 1 

41 


100-150 


3 


2.6342(56) 




51 


100-150 


3 


2.6741(32) 




61 


inn 1 

100-150 


o 

3 


2.6977(43) 




71 


1 AA 1 P' O 

100-150 


3 


O TO O z"" 1 / OO \ 

2.7036(39) 




101 


100-150 


3 


2.7234(25) 




oo 






2.74(1) 


0.2275 


21 


50-250 


8 


2.8858(56) 




31 


50-150 


5 


2.9938(40) 




41 


~i o o 1 rf\ 

100-150 


3 


3.0365(51) 




51 


100-150 


3 


3.0578(37) 




61 


100-150 


3 


3.0725(49) 




71 


100-150 


3 


3.0832(32) 




oo 






3.10(1) 


0.23 


11 


16-30 


4 


2.7649(89) 




15 


16-46 


7 


O 0,1 ,1 O / O /"I \ 

3.0443(30) 




O 1 

21 


OO A OO 

30-400 


1 o 

12 


O 1 TO/^/ooN 

3.1726(20) 




25 


30-80 


6 


3.2200(28) 




31 


30-90 


7 


3.2608(17) 




35 


40-120 


5 


3.2790(35) 




41 


40-120 


5 


3.2904(22) 




51 


76-200 


6 


3.3034(13) 




61 


60-150 


4 


3.3118(14) 




81 


40-120 


5 


3.32151(95) 




101 


50-100 


6 


3.32622(58) 




119 


50-80 


4 


3.32875(66) 




151 


50-90 


5 


3.33083(40) 




171 


50-90 


3 


3.33128(57) 




oo 






3.335(3) 
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K 


L 


H 


N H 




0.235 


31 


50-90 


5 


3.5859(17) 




35 


50-90 


5 


3.5939(11) 




41 


50-90 


5 


3.6029(18) 




45 


50-90 


5 


3.6073(13) 




51 


rn f\f\ 

50-90 


5 


3.6132(14) 




61 


50-90 


5 


3.62088(83) 




81 


50-90 


5 


3.62346(67) 




oo 






3.630(5) 


0.24 


11 


16-30 


4 


3.5778(22) 




15 


16-46 


7 


3.6852(34) 




21 


20-400 


13 


3.7447(14) 




25 


30-80 


6 


3.7616(26) 




31 


30-90 


7 


3.7829(19) 




35 


20-120 


6 


3.7875(15) 




41 


20-120 


6 


3.7943(13) 




51 


26-200 


8 


3.79888(24) 




61 


30-210 


7 


3.80299(95) 




81 


40-120 


5 


3.80849(81) 




101 


50-100 


6 


3.81055(57) 




151 


50-90 


5 


3.81181(32) 




oo 






3.813(2) 



30 



K 


L 


H 


N H 




0.245 


31 


50-90 


5 


3.9101(16) 




35 


50-90 


5 


3.9141(13) 




41 


50-90 


5 


3.9181(10) 




45 


50-90 


5 


3.9202(16) 




51 


50-90 


5 


3.92080(67) 




61 


50-90 


5 


3.92691(61) 




81 


50-90 


5 


3.92687(63) 




oo 






3.930(2) 


0.25 


11 


16-30 


4 


3.8552(26) 




15 


16-46 


7 


3.9297(21) 




21 


20-400 


13 


3.96912(78) 




25 


20-80 


7 


3.9839(16) 




31 


30-90 


7 


3.9906(15) 




35 


20-120 


6 


3.9929(12) 




41 


20-120 


6 


4.0028(11) 




51 


26-200 


8 


4.00478(52) 




61 


30-210 


7 


4.00704(83) 




81 


40-120 


5 


4.00900(50) 




101 


50-100 


6 


4.01017(32) 




151 


50-90 


5 


4.01138(20) 




oo 






4.012(1) 



Table 6: The estimated values of e s for some values of K and L are shown. 
The column H shows the height region over which the average is calculated. 
Nh is the number of different H used in the averaging operation. The values 
of e s in L = oo denotes the estrapolated values of e s (K,L) to L = oo and 
they are the estimates for the interface energy in the thermodynamics limit. 
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year a 

1982 1.05 ±0.05 
1985 1.2 ±0.1 

1988 1.58 ±0.05 
1993 1.5 ±0.1 

1993 1.25 ~ 1.68 

1993 1.52 ±0.05 

1993 1.221 ~ 1.492 

1993 1.92 ±0.15 



reference 



5 

II 
13 



m 



1993 1.42 ± 0.04 present study 



Table 7: The estimates of o"o are shown with the year of publication. 



n 


A(n),256 3 


A(n), 1024 3 


2 


2.4625(1) 


2.4624(2) 


3 


2.4738(5) 


2.4735(5) 


4 


2.481(1) 


2.480(2) 


5 


2.485(3) 


2.485(5) 



Table 8: The estimates of A(n) are shown. The estimates from 256 3 and 
1024 3 lattice coincide within the present accuracy of the estimation. 
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